c
c    To be linked with wp_xi_cf_hod_12.f
c    This module gives the N1(M),N2(M) and 1-halo term power spectrum
c    designed to model the AGN N1(M) 
c    N1(M) is usual Truncated power-law and assumed to be all satellites.
c    
c    N2(M) is taken from an ascii table gives separate entries for central 
c       and satellite galaxies.
c      
      subroutine init_nm(argc2)      
      ! Read model parameters/ In this model
      ! Mcrit_1,alpha1,filename
      parameter (nmass=50,idv=70)
      character*(*) argc2
      character c_or_s*1,nm2_fil*80,rline*132,cut_space*132      
      integer nm2(0:1)  ! index=1 centrel, 2: satellite
      real lgm2(nmass,0:1),gn2(nmass,0:1),lgmx,gnx
      common/wgt2par/nm2,lgm2,gn2
      common/wgt1par/rm_crit1,alpha1,f_cen1
      !
      read (argc2,*) rm_crit1,alpha1,f_cen1,nm2_fil
      ! For input, if rm_crit=0 means no such population. In the linear
      ! scale, it is converted to 1.0 so that its log will give 0.
      if (rm_crit1.gt.0.and.rm_crit1.lt.1e5) then
         rm_crit1=10**rm_crit1
      else if (rm_crit1.le.0) then
         rm_crit1=1.0
      end if
      !
c
c     Read Table for population 2
c
      open (unit=70,file=nm2_fil,status='old')
      nmc = 0
      nms = 0
      do while (.true.)
         read (70,'(A)',end=3) rline
         rline = cut_space(rline)
         if (index(rline,'!').ne.1) then
            read (rline,*) c_or_s,lgmx,gnx
            if (lgmx.gt.1e5) lgmx=log10(lgmx) !If entry is lin, covert to log
            if (c_or_s.eq.'c') then
               icen=1 ! Central
            else
               icen=0 ! Satellite
            end if
            nm2(icen)=nm2(icen)+1
            lgm2(nm2(icen),icen)=lgmx
            gn2(nm2(icen),icen)=gnx
         end if
      end do        
 3    close (70)
      ! Output Header
      write(6,'(A)')
     $     '! CCF model pl_tbl/Sample1: Truncated PL w/ central'//
     $     ' frac=f_cen1 Sample2: Table'
      write(6,'(A,2(1X,F6.2),F7.3)')
     $     '! Sample1: Log(M_cr1) alpha1 f_cen1:',
     $     log10(rm_crit1),alpha1,f_cen1
      write(6,'(A,A,1X,I2,A,I2,A)')
     $     '! Sample2: N(M) file:',nm2_fil(1:index(nm2_fil,' ')-1),
     $     nm2(0),' rows (sat) ',nm2(1),'rows (cen)'
      end
c
c
      real function galn_halo (rm,isam,icen) 
      ! rm: Mass of DMH in solar mass
      ! isam: sample number
      ! icen=0: satellite galaxy, icen=1: central galaxy 
      parameter (nmass=50)
      integer nm(0:1)
      real lgm(nmass,0:1),gn(nmass,0:1),lgmx
      common/wgt1par/rm_crit1,alpha1,f_cen1
      common/wgt2par/nm,lgm,gn
      ! By definition, if rm_crit*<=10, no such population.
      ! This controls whether it is central only or satellites
      ! only... But... this is not quite right...
      ! center have to be l<=1. Need correction if central object
      ! is considered. Currently all AGNs are assumed to be satellites.
      if (isam.eq.1.and.icen.eq.0) then ! Sample 1 Satellite
         if (rm.lt.rm_crit1.or.rm_crit1.le.10) then             
            galn_halo = 0
         else
            galn_halo = ((rm/rm_crit1)**alpha1)*(1.0-f_cen1)
         end if
      else if (isam.eq.1.and.icen.eq.1) then ! Sample 1 Center
         if (rm.lt.rm_crit1.or.rm_crit1.le.10) then 
            galn_halo = 0
         else
            galn_halo = ((rm/rm_crit1)**alpha1)*f_cen1
         end if
      else ! Sample 2
         lgmx = log10(rm)
         if (lgmx.lt.lgm(1,icen)) then
            galn_halo = 0
         else
            imin=1
            imax=nm(icen)
            if (lgmx.gt.lgm(nm(icen),icen)) then
               imin = imax-1
               ff=1.0
            else
               call intra_mono(lgm(1,icen),lgmx,imin,imax,ff)
            end if
!      
            galn_halo =  (1.-ff)*gn(imin,icen)+ff*gn(imax,icen)
            if (icen.eq.0) galn_halo = max(galn_halo,0.) ! Sate
            if (icen.eq.1) galn_halo = min(max(0.,galn_halo),1.0) ! cen        
         end if
      end if      
      return
      end
c
c
      real function wgt_cf(rm,icen) ! Calculate Weight for
                              ! Expected number of pairs 
                              ! between sample 1 and sample 2
      rm10 = rm/1e10          ! rm in units of 1e10 solar mass/to match with wgt_samp   
      if (icen.eq.0) then     ! Sat vs Sat
         wgt_cf = galn_halo(rm,1,0)*galn_halo(rm,2,0)/(rm10*rm10)
      else   ! Central vs satellite
         wgt_cf = (galn_halo(rm,1,1)*galn_halo(rm,2,0)+
     $        galn_halo(rm,1,0)*galn_halo(rm,2,1))/(rm10*rm10)
      end if   
      return
      end
c
      subroutine model_help()
      print '(A)',' Model Specification: ccf_pl_tbl2_cut.f'
      print '(A)',' N(M)_pars: Mcr1,alpha1,f_cen1,file2'
      print '(A)','  Sample1: N1(Mh) =0 (Mh<Mcr1_s,Mh>Mcut); '//
     $     '\\propto(Mh/Mcr1)^alpha1 (Mcr1_s<Mh<Mcut) (sat+cen)'
      print '(A)','  Sample1: N1_c(Mh)=f_cen1*N1(Mh)  (cen)'
      print '(A)','    Mcr=0 means zero N(M) for the population.'
      print '(A)','    Can give Log10(Mcr) instead of Mcr (automatic).'
      print '(A)','  Sample2: ASCII file ''file2''  with 3 cols:'
      print '(A)', '  LogMH [h^{-1}M_sun],N(M)_satellite,N(M)_center'
      print '(A)', ''
      return
      end
c
